import pandas as pd
import numpy as np
from scipy.stats import norm
import statsmodels.api as sm
import warnings

def topsis_ranking(df: pd.DataFrame, weight_series: pd.Series, output_path: str, num_levels: int = 4) -> tuple:
    """
    使用TOPSIS方法进行排名和分档分析。
    
    参数:
    df: 包含评价指标的DataFrame，必须包含权重计算所需的所有列
    weight_series: 各指标的权重，索引对应的列名必须与df中的列名一致
    output_path: 结果保存的CSV文件路径
    num_levels: 分档等级数量
    
    返回:
    tuple: 包含处理后的DataFrame，分界点列表和模型摘要文本
    """
    # 验证输入类型
    if not isinstance(df, pd.DataFrame):
        raise TypeError("df必须是pandas DataFrame类型")
    if not isinstance(weight_series, pd.Series):
        raise TypeError("weight_series必须是pandas Series类型")
    if not isinstance(output_path, str):
        raise TypeError("output_path必须是字符串类型")
    
    # 检查必须存在的列
    if not all(col in df.columns for col in weight_series.index):
        missing_cols = [col for col in weight_series.index if col not in df.columns]
        raise ValueError(f"df中缺少必要的列: {missing_cols}")
    
    # 数据预处理和Topsis算法实现
    try:
        indicators = weight_series.index.tolist()
        matrix = df[indicators].values
        # 使用更安全的数值计算方法
        norm_matrix = matrix / np.sqrt((matrix ** 2).sum(axis=0)[:, np.newaxis])
        
        # 使用clip避免极端值影响
        weight_values = np.clip(weight_series.values, 1e-10, 1.0 - 1e-10)
        weighted_matrix = norm_matrix * weight_values
        
        ideal = np.max(weighted_matrix, axis=0)
        anti_ideal = np.min(weighted_matrix, axis=0)
        
        dist_pos = np.sqrt(((weighted_matrix - ideal) ** 2).sum(axis=1))
        dist_neg = np.sqrt(((weighted_matrix - anti_ideal) ** 2).sum(axis=1))
        
        # 处理可能的除零情况
        with np.errstate(divide='ignore', invalid='ignore'):
            closeness = dist_neg / (dist_pos + dist_neg)
        np.seterr(divide='ignore', invalid='ignore')
        closeness = np.where(np.isnan(closeness), 0.5, closeness)
        
        df['相对接近度C'] = closeness
        
        if '年份' in df.columns:
            df['年度排名'] = df.groupby('年份')['相对接近度C'].rank(method='min')
            
    except Exception as e:
        raise RuntimeError(f"TOPSIS计算过程中出错: {str(e)}")
    
    # Probit分档分析
    try:
        epsilon = 1e-10
        df['C值'] = df['相对接近度C'].clip(epsilon, 1 - epsilon)
        
        # 当有效样本量较小时，使用Probit变换
        if len(df) < 5:
            df['C值'].replace([0, 1], [epsilon, 1-epsilon], inplace=True)
            
        df['Probit值'] = norm.ppf(df['C值'])
        
        # 为OLS回归准备数据
        X = sm.add_constant(df['Probit值'])
        y = df['C值']
        
        # 添加小噪声以防止完全拟合，提高泛化能力
        y_noisy = y + np.random.normal(0, 1e-4 * (y.max() - y.min()), len(y))
        
        model = sm.OLS(y_noisy, X).fit()
        df['C拟合值'] = model.predict(X)
        
    except Exception as e:
        raise RuntimeError(f"Probit分析过程中出错: {str(e)}")
    
    # 根据分档级别设置概率阈值
    if num_levels == 4:
        probs = [0.25, 0.5, 0.75]
    elif num_levels == 3:
        probs = [1/3, 2/3]
    elif num_levels == 5:
        probs = [0.2, 0.4, 0.6, 0.8]
    else:
        raise ValueError("不支持的分档级别，只支持3、4、5级分档")
    
    # 使用模型预测确定分界点
    try:
        probit_thresholds = [norm.ppf(p) for p in probs]
        const, slope = model.params
        
        # 计算分界点时使用有符号规范化
        cut_points = const + slope * np.array(probit_thresholds)
    except Exception as e:
        raise RuntimeError(f"计算分界点时出错: {str(e)}")
    
    # 分档处理
    try:
        # 创建分界点列表
        bins = [-np.inf] + list(cut_points) + [np.inf]
        labels = [str(i) for i in range(1, num_levels + 1)]
        
        # 使用pd.cut进行分档
        df['分档等级'] = pd.cut(df['C拟合值'], bins=bins, labels=labels, right=True)
        
    except Exception as e:
        raise RuntimeError(f"分档过程中出错: {str(e)}")
    
    # 输出结果
    try:
        df.to_csv(output_path, index=False, encoding='utf-8-sig')
    except Exception as e:
        raise RuntimeError(f"保存输出文件时出错: {str(e)}")
    
    return df, cut_points.tolist(), model.summary().as_text()
